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Abstract 

We present a new algorithm, the Time Dependent Phase Space Filter 
(TDPSF) which is used to solve time dependent Nonlinear Schrodinger 
Equations (NLS). The algorithm consists of solving the NLS on a box 
with periodic boundary conditions (by any algorithm). Periodically in 
time we decompose the solution into a family of coherent states. Coherent 
states which are outgoing are deleted, while those which are not are kept, 
reducing the problem of reflected (wrapped) waves. Numerical results are 
given, and rigorous error estimates are described. 

The TDPSF is compatible with spectral methods for solving the in- 
terior problem. The TDPSF also fails gracefully, in the sense that the 
algorithm notifies the user when the result is incorrect. We are aware of 
no other method with this capability. 



1 Introduction and Definitions 

Consider a semilinear Schrodinger equation on M^+-'^ 

idtibix, t) = -{l/2)^i,{x, t) + g{t, X, ^{x, t))^{x, t) (1.1) 

where g{t,x,-) is some semilinear perturbation. For instance, g{t,x,-) could 
be V{x,t) + f{\tlj{x,t)\^) for some smooth function / and (spatiaUy) locahzed 
potential V{x,t). Abusing terminology, we refer to g{t,x, •) as the nonlinearity 
even when it is linear. 

We assume the initial condition and nonlinearity are such that the nonlin- 
earity remains localized inside some box [— Lnl, ^nl]^- Outside this region the 
solution is assumed to behave like a free wave. 

One very common method of solving such a problem is domain truncation. 
That is, one solves the PDE p.ip numerically on a region [— L, L]^. On the 
finite domain, boundary conditions must be specified. Dirichlet or Neumann 
boundaries introduce spurious reflections, while periodic boundaries allow out- 
going waves to wrap around the computational domain. This causes the numer- 
ical solution to become incorrect after a time T « L/fcmax, where fcmax is the 
"maximal velocity" (the highest relevant spatial frequency) of the solution. 



1 



There are two major approaches to deahng with this problem. The Dirichlet- 
to-Neumann method consists of attempting to use the exact solution as a bound- 
ary condition. This works rather well for the wave equation [J [51 [51 [TTl [TB] . 
but it is fraught with problems for dispersive waves and only limited progress 
has been made [111 [IHl UHl [2S] ■ Dirichlet-to-Neumann methods also prevent the 
use of fast spectral methods for solving the interior problem because spectral 
methods based on the FFT naturally impose periodic boundary conditions. A 
spectral solver for the interior problem is highly desirable since dispersive waves 
are difficult to calculate by other methods (see Remark [XT]) . 

The other approach is to add a dissipative term which is localized on a buffer 
region [— (Ljnt + w), (Lint + w)]^ \ [—Lint, ^int]''^- The dissipative term can be 
an absorbing potential [17] or the more sophisticated PML [3] . This dissipates 
outgoing waves, but also dissipates incoming waves located near the boundary. 
This can introduce spurious dissipation (see Section [372]) which does not occur 
for the TDPSF. 



1.1 Our Approach 

We propose a new method. Suppose we want to approximate the solution in 
the space (with s large enough to make (jl.ip is well posed). We make the 
assumption that outside some region [—Lp^Lp]^, ip{x,t) behaves like a free 
wave, I.e. size of the computational box Ljnt 

must be larger than Lp. 

We assume the nonlinearity is nearly localized both in position and momen- 
tum. That is, we assume the existence of Xnl, fcmax,NL so that: 

\\[^- X[-L^l^.L^l^]«{x)]g{t,x,^P{x,t))'^P{x,t)\\^^ « (1.2a) 

[1 - X[-fc„,,,NL,fc„ax,NL]"(fc)]5(i,^X2;,*))V'(a;,0 „ ~ (1.2b) 

For the rigorous proof, we assume that the nonlinearity is Lipschitz in H'^, 
though this is probably unnecessary. In particular, the nonlinearity could take 
the form of a time dependent short range potential V{x,t)'ip{x,t) or a local 
nonlinearity f{\'ip{x,t)\ '')'ip{x,t) with f{z) a bounded function. 

We also assume that the solution remains localized in frequency, that is the 

norm of V'(fc, t) is small outside [— fcmax, fcmax]^ for some large number fcmax 
(the maximal momentum of the problem, which we assume exists) . 

Our algorithm is as follows. We assume the initial data is localized on a 
region [-Lint, Lint] ^. We solve (|1.1|) on the box [-(Lint + w),Lint -I- w]^ on 
the time interval [0, Tstcp]- We insist that Lint > Lf,L^i^ so that the solution 
behaves like a free wave on [—(Lint + w), (Lint + w)]^ \ [— Li„t, Lint]^- 

By making Tstcp smaller than w/Sfcmaxj we can ensure that -0(2^, LLtcp) is 
mostly localized inside box [—(Lint -I- w), (Lint + w)]^. Since very little mass has 
reached the boundary, the error is nearly zero [231 Thm. 4.1]. 

We then decompose il^^^, Lgtep) into a sum of Gaussians (indexed by a, 6 G 
Z^, with lattice spacing xq in position, in momentum and standard deviation 
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a): 



\aXo\ao<Liat + W 
|&fco|oo<fcmax 

We then examine the Gaussians near the boundary, |aa;o|oo > ^int, and deter- 
mine whether they are leaving the box or not under the free flow {\x\p denotes 
the IP norm on finite vectors). By assumption, e'^^/^-''^* is a sufficiently accurate 
approximation to the true solution for this part of the solution. We apply this 
approximation to each framelet: 

^i(l/2)At^-N/4^-N/2^ikob-S^-\S~axo\l/2cr^ 

^ cxp (ibko ■ {x-{bko/2)t - qxq)) f -\x-bkot-dxo\l \ 

~ 7r^^/%^/2(l + ii/a2)A'/2 ^"^^ y 2a^{l + it/cr^) J ^ ' 

Under the free flow, a framelet moves along the trajectory axQ + bkot while 
spreading about it's center. If a given Gaussian is leaving the box, we delete it 
(set ip^^ =0). Some Gaussians spread more quickly than their center of mass 
moves, and we do not present here an algorithm to deal with these Gaussians. 

After this filtering operation, the only Gaussians remaining are either inside 
the box [—-Lint, -t'int]^, or inside the box [—(Lint + w), (iint + w)]^ but moving 
towards [—Lint, Lint]^ ■ Thus, it is safe to propagate for a time Tstep, since what 
remains will not hit the boundaries before this time. We do the same at time 
2Tstcpj STstopj etc. 

The main drawback to this approach is that some Gaussians are ambiguous. 
Consider a Gaussian with velocity 0: Tr^^/^cr^^/^g-lai-aajols/ao- _ Suppose also 
that axo is located inside [— (iint+w), (Lint+w)]^\[— iint, Lint]^- By examining 
(|1.3p . we observe that this framelet is spreading out laterally both into the box 
and outward. If we delete it, we have removed waves which should have returned. 
If we fail to delete it, then that part of the framelet which is spreading outwards 
will wrap around and cause an error. 

This happens only for Gaussians which are moving slowly relative to the 
box. Thus, we impose one additional assumption - we assume that Gaussians 
moving so slowly that they move both ways do not occur in the solution. The 
algorithm to deal with low frequency waves is more involved and we relegate 
this to a future work |23j . 

1.2 Error Bounds 

We prove rigorous error bounds for the TDPSF algorithm in [23], as well as 
stating explicitly the assumptions and defining what we mean by . Although 
the error bound is too long to state here, we describe briefiy the form it takes. 
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For a general time-stepping algorithm (with periodic boundaries and no 
filtering), the error bound would take the following form: 

sup \\U(t)iJo{x) - *(a;, < BoundaryError(rmax) 

te[o,T„ax] 

+ HighFrequency(Ti„ax) + LowFrequency(ri„ax) 

+ NonlocalNonhnearity(r,nax) + Instabihty(Tmax) (1-4) 

The term BoundaryError(Tmax) encompasses errors due to waves wrap- 
ping/reflecting from the boundaries of the box. For many problems, this is the 
dominant error term. It is directly proportional to the mass which would have 
(if we were solving the problem on M^) radiated outside the box [—Lint, ^int]''^- 

The HighFrequency(rmax) part stems from waves with momenta too high to 
be resolved by the discretization. The term LowFrequency(Tniax) encompasses 
errors due to waves with wavelength that is long in comparison to the box. The 
term NonlocalNonlinearity(rmax) stems from that fraction of the nonlinearity 
itself which is located outside the box. The Instability(rmax) stems from the 
possibility that the dynamics of the solution itself might amplify the other errors 
dramatically (e.g. in strongly nonlinear problems). 

Our algorithm reduces the term Boundary Error(Tniax)- The other errors are 
assumed to be small, since we are concerned only with the boundary error. 

The method of proof is a direct calculation. We calculate the errors made 
in the filtering step as well as the errors made while propagating in between 
filtering. We then add up the error over time, taking into account instabilities 
of the system being simulated. The error takes the form described in (|1.4[) . 
The Boundary Error (Tniax) term can be reduced by increasing w, and the error 
behaves like BoundaryError(Ti„ax) — 0(e~'^"'T„iax)- In some cases (problems 
which are asymptotically complete and have no bound states), we believe that 
this error can be proved to be time independent. 

The main drawback of our algorithm is that it does not provide us the ability 
to filter waves for which the wavelength is longer than the buffer region. This is a 
due to the Heisenberg uncertainty principle. This problem will be addressed by 
a novel multiscale argument in [23| . The problem is shared with most methods 
of absorbing boundaries. The TDPSF algorithm also alerts the user when it 
fails, unlike all other methods we are aware of. 

2 The Algorithm 

2.1 The Windowed Fourier Transform 

The discrete windowed Fourier transform frame is a way of localizing a function 
f{x) in the phase space. For the big picture, see [6l[7l[8]; see also [24t Section 
3] for technical details needed by the TDPSF. 

The WFT is an expansion of a function f{x) into the following "basis" 
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(actually a frame) for L'^{R'^): 

^ia.b)i^) = 7r-^/V"~/2e''=«''-^e-l^-'^"«l^/2"', a,beZ^ x (2.1) 

The lattice spacing in position and momentum, xq, ko, and the standard devia- 
tion a must positive real numbers. To be a frame, it must hold that xok^ < lix. 

This family of functions is not a basis in the usual sense, but a frame. A 
family of functions {0j(a;)}jgj is a frame in the Hilbert space H if there exists 
an < Ap < Bp < cxd so that: 

Af \\f{x)\\H < (^KflM^))]'^ < Bp \\f{x)\\^ 

A frame is an overcomplete, non-orthonormal basis with a numerically stable 
reconstruction procedure. We can then decompose any / e X^(R^) as: 

/(^) = E /(H,?)<^(a-i)(^) (2-2) 

fia,t) = (fi^MiaM^)) (2-3) 

The functions i'l^^l^ix) are given by 0^-g^(a;) = e^'^°^'^g{x — axo) with g{x) £ 
L^(R^) (see [SlIZ] for procedures to calculate g{x)). 

In [24j we extend the analysis of [7j and show that if xofco = Stt/M for 
M £ 2N, g{x) decays exponentially in \x\i and g{k) decays exponentially in 
\k\i: 

\d:m\ <g(.^o,fco,A^,a)e-"-(^°^^°)l"^l^ (2.4) 

g{xo, ko, N, a) is a constant (explicitly bounded in [24 ) as is r(a;o, fco) = xqM/Sttct. 
(|2.4p holds for g{k), but with xq and fco swapped and replacing a. 



2.1.1 Phase Space Localization 

The WFT allow us to define a concrete realization of phase space. From here 
onward, we will consider x to be a discrete realization of phase space. 
The vector (a, b) G Z^ x will represent the point at axQ in position, and bko 
in momentum. 

With this in mind, we can now construct phase space localization operators 
very simply. For a set F e x Z^, we define Vp by: 

Vpi,{x)= ^ = E {^iarb)(^)\^(^)) 'l^iaM^) (2-5) 

(o,b)e_F {a,b)eF 

It can be shown [SJ [B] that phase space localization in terms of the WFT 
corresponds closely (though not exactly) to the usual phase space localization. 
We state one theorem (proved in [24], based on a result by Daubechies from 
[3 inj) as an example of this: 
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Theorem 2.1 Let Bx = [-X,X]^, Bk = [~K,K]^ for X,K <oo. Suppose 
also that xoko = 27r/M for M G 2N. Then there exists C , X and K so that if 
X' ^ X ~X, K' ^ K -K, then: 

\\f{x) -''PB^,xB,,,f{x)\\^, 

< C (||(1 - P|,;,„(x))/(f)||^, + ||(1 - + 6 11/11^.) (2.6) 

The constant C depends on e, s, a;o, ko, a, X and K , as do X and K. C is given 
explicitly in 



2.1.2 Computation of the WFT Coefficients and Phase Space Pro- 
jections: How to do it, and how hard it is 

We now present the algorithm for computing the framelet coefRcients. The 
algorithm consists of computing the products f{x)g{x—ax{)), followed by Fourier 
transforming the results. Due to the spatial decay of g{x), we can truncate the 
domain to the box [— Le, Le]^ with minimal error. 

Note that g(x) can be calculated as accurately as necessary [SlIZ], and we 
will not discuss this here. 

Algorithm 2.1 Calculation of Windowed Fourier Transforms 

This algorithm calculates, for a function f{x), the WFT coefficients fj^^ for 

a € A C . We assume that the frequencies are hounded above by /cmax • The 
lattice spacing in frequency, ko is taken to be 2i:/ L^; taking it to he any larger 
than this yields only logarithmic improvements in computational complexity. 

1. Let AcZ^ be some set of position coordinates. 

2. For each a ^ A multiply f{x)g{x — axo) for x £ [—L^, L^]^ + oxq only. 

3. Calculate the Fast Fourier Transform of f{x)g{x — axo) on this region. 
For each a, the resulting function is f^^ . 

We observe that this algorithm is local in space. This means that if A is 
finite, then the computational cost is proportional to: 

\A\ kZ,. |supp5(f)r log(fc„,ax Isuppg(f)l) = \A\ 2^fc,l,Lf 

The reason for this complexity is as follows. For each a S A, we need to 
compute an FFT. The FFT is computed only on the "support" (after truncating) 
of g{x) and the lattice spacing in this region is 27r/fc,nax = 0(l/fcniax)- Thus, 
there are 0(|supp5(x)| fc^ax) = -^^ data points in this region. The FFT has 
computational complexity Mlog(M). In addition, we need to compute \A\ of 
these FFT's. 

Note that this algorithm can be parallelized very easily, simply by having 
different processors compute the FFT's for different ax^. 
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Remark 2.2 As an example, consider the case when A = {a Q : axo G 
[—(Lint + w), (Lint + w)]^ \ [—Lint, ^int]''^}- Then the size of A is proportional 
to I [-(Lint + w), (Lint + w)]^ \ [-Lint, Lint]^ | / . If Lint > w, then this is of 
order L^^~^, and the computational complexity is 0{L^^^ k^^^\og{Lintkma.x)) ■ 
We consider this case since this is what the TDPSF requires. 

Phase space projections can also be computed. Let F C x be finite. 
Let A = {a € : 3b £ Z,{a,b) G F}. Then we provide the phase space 
projection algorithm: 

Algorithm 2.2 Phase Space Projection Algorithm 

This algorithm computes the phase space projection onto a region of phase 
space F. 

1. Compute A = {a : (a, h) G Assume this is finite. 

2. Compute /(^gj for a £ A as in Alaorithm \2.1\ 

3. Define a new function /* : Z^ x Z^ C hy: 

^. Compute the inverse WFT of f^ . The result approximates Vpfix), with 
errors caused by the truncation of g{x) in Alaorithm \2.1[ 

Clearly the computational complexity of Algorithm 12.21 is of the same order 
as that of Algorithm 12. II 

2.2 Propagation with Periodic Boundaries: FFT/Split Step 
Algorithm 

The TDPSF is a filtering algorithm, which is built on top of another propagation 
method which solves on a box with periodic boundaries. The exact manner 
in which this is done is irrelevant for our purposes, provided it is sufficiently 
accurate. However, it is important to note that the spectral propagation method 
is significantly better than the usual FDTD methods, and that compatability 
with spectral methods is an important feature of the TDPSF. 

Fix a grid spacing Sx and timestep 6t. The computational grid is sim- 
ply the region [— Lcomp, Lcomp]^ sampled at the lattice with spacing Sx. This 
corresponds to a lattice spacing in momentum of 27r/Lcomp, with maximal mo- 
mentum 2Tr/6x. A common rule of thumb is that if the problem has a maximal 
momentum fc,nax, then 6x = 47r/fcniax- 

Algorithm 2.3 Split Step Propagation Algorithm 
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This algorithm approximates the propagation of (jl.ip on the region 
[— Lcomp, -Z^comp]^- The accuracy is 0{5t^) time and 0{5x^) ("spectral accu- 
racy") in space. It is 0{St^) in time if (jl.ip is linear. 

First, the discrete operator e*'^/^^'^*/(a;) is defined by computing the FFT of 
f{x), multiplying by e'^^ * and then computing the inverse FFT of the result. 

With this in mind, the algorithm is as follows (taking ipo{x) as initial con- 
dition, and assuming t is a multiple of 6t): 

1. Define V'i/2(a;) = e'(i/2)A5t/2^Q(a;). 

2. For j ^0,..., t/St - 2, define: 

3. Finally, define ij^t/stix) ~ il}{x,t) by: 

This is a numerical realization of the Trotter-Kato product formula: 



Uit) 



,,(l/2)A(-5t/2) 



't/St 



/2)A5t^~ig{j5t,x,ip{x,jSt))6t 



J{l/2)ASt/2 



(2.7) 



The computational complexity of Algorithm 12.31 is 

0((Lint + W)^fc^axln[(iint + w)kmax]) 

per timestep (since the number of data points is of order (Ljnt + w)fcniax)- The 
number of timesteps is T^^^/ 5t. Thus, split step propagation has time com- 
plexity 0[I/^t^^axlog(ii„tfci„ax)(T,„ax/(5i)]. Algorithm [13] is by now a textbook 
result 0], so we will not discuss it further. 



2.3 The TDPSF Algorithm 

The TDPSF algorithm consists of solving (jl.ip on the box [— icomp, -^comp]^ 
with periodic boundary conditions, using Algorithm 12.31 or any other method. 
At times Tgtep, 2Tstep, etc, we subtract from the solution the projection onto 
Gaussians which are outgoing by means of Algorithm 12.21 

Outgoing Gaussians are defined to be those (a;) with oxq G [— (^int + 
w), (Lint + w)]^ \ [— iint, ^int]^ and which are moving strictly outward under 
the free flow; that is, those which are leaving [—Lint, ^int]^ soon (before a time 
Tstcp has passed) and which will never return. 

Incoming Gaussians are those which will not leave [— i^comp, ^comp]^ before 
time Tstep has passed. 

There is also a set of ambiguous gaussians. These are gaussians with ve- 
locity sufficiently low so that they spread about their center faster than they 
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move. They are ambiguous because they are both leaving [— Lcomp, Lcomp]^ 
and returning to [— Lint, ^int]''^ due to their spread. An example of this would 
be e~(^"^'"'~"'/^) Z^"' . This gaussian is located inside the buffer region, and is 
spreading laterally with no motion whatsoever. 

The TDPSF algorithm does not know what to do with these gaussians, 
which correspond to low frequency waves. However, unlike other algorithms (e.g. 
Dirichlet-Neumann, PML or absorbing potential), the TDPSF does know when 
such waves are causing a problem. Thus, incorrect results are never returned. 

Algorithm 2.4 Propagation algorithm This is the main propagation algorithm. 
In this algorithm, we consider a, xq, ko to be fixed. Lint cind T^ax cife also con- 
sidered fixed. The initial data is considered fixed, and localized inside [—Lint, Lint]^ 
The approximation will be denoted by 'i!(x,t). Also, fix a small tolerance e > 0. 

1. Before beginning, precalculate the set of framelets which are outgoing, and 
those which are ambiguous. 

2. Define ^(a;,<) iteratively as follows. Loop over n = . . . Tmax/Lstop- In 
what follows, the propagator U(t) is calculated by Alaorithm \2.3[ 

(a) Fort G [{n - l)Tstcp, "Tstop), define 

ip{x,t) ^U{t- {n- l)Tstcp)'4'{x, (n - l)Tstcp) 

(b) For t — nTstcp, define: 

"^{x, nTstcp) = (1 - 'PovT)l^{TstcpH{x, (n - l)Tstep) 

(c) At each integer multiple o/Lgtcp, compute 

\\VMABU{Tstcp)lp{x, (n - l)'rstcp)|l^s 

where AMB is the set of ambiguous framelets. If this quantity is 
greater than e, stop the program and notify the user. 

Remark 2.3 The TDPSF fails gracefully in the following respect. If the TDPSF 
is unable to filter due to slowly moving waves, the program fails and the user is 
notified. Absorbing potentials, the PML and Dirichlet-Neumann maps do not 
have this property. 

Remark 2.4 In order to calculate "Pout, we must calculate the WFT coeffi- 
cients of the solution in the region [— Lcomp, Lcomp]^ \ [—Lint, Lint] ^. To cal- 
culate Pamb, we need the WFT coefficients from the same region. Thus, only 
one WFT is needed for steps (b) and (c). The cost of computing the norm 
of the ambiguous framelets is cheaper than the cost of a WFT; this is merely a 
sum rather than a set of FFT's. 
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Thus, all that remains is to do is explain which framelets are outgoing, 
incoming and ambiguous. 

For a given (a, 6), this can be determined easily, by a direct calculation based 
on (|1.3p . This is done in for s = 0, 1, and other cases can be done with the 
help of Maple/Mathematica. 

The end result of the calculation is the following. Define by 



if we measure the error in i^f 



(2.8a) 



R^{t) ^ y/a^+t^/a^ max ■ 



N/2, 



e^7r^/2 



2|5W-i|(l 
{N + 2)/2 



Ifefcoli), 



215^- 



-,1/2 



1/2 



(2.8b) 



if we measure the error in H^{M.^). 

The function T~^{a,x) is the inverse function of a; i— > r(a,x), with T{a,x) 
the complementary incomplete Gamma function. This definition implies that 
all the mass (up to a tolerance e) of e''^^/^^'^*(/)jg (x) is, at time t, contained in 

a ball of radius i?g about the point axo + bkot. That is: 



^i{l/2)At^-N/4^~N/2^ikob-x^-\x-ax„\l/2a 



_H"»i(-B(t)^) 



< e 



B{t) = {x : \axo + bkot - x\2 < R^{t)} 



(2.9a) 



(2.9b) 



Thus, if B{t) does not intersect [—Lint, Lint] , this framelet is strictly bad. 

Similar calculations can be done in if" , although i?g must be redefined (and 
now depends on b). Thus, given the result of this calculation, we define the set 
of outgoing Gaussian's to be the set of (a, b) G x such that: 



1. aXQ e [-(Lint + w), {Lint + W)]^ \ [-Lint, Lint]^ 

2. B{t) does not intersect [— Lint, Lint]^ for any t < Tmax- 

A simple criterion to determine whether a given framelet is outgoing is to 
check whether axg G [—(Lint + w), {Lint + w)]^ \ [— Lint, Lint]^, and whether 
bnko > 2 In e/cr. Here, 5„ is the component of bko pointing in the direction nor- 
mal to [—Lint, Lint]^. Conversely, this means that to filter waves with frequency 
fcmin (or higher) making an error e, then a > Ine/fcniin- This determines the 
width of the buffer region. This appears to be somewhat larger than what a 
PML requires, although we were unable to obtain high accuracy using current 
implementations of the PML (see Remark l3.ip . 

The set of ambiguous framelets is the set of framelets for which B{t) inter- 
sects [— Lint, Lint]^, but also intersects \ [—(Lint + w), {Lint + w)]^ . That is 
to say, these are framelets with both incoming and outgoing components. 
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A calculation based on (|1.3p shows that these consist of framelets with either 
low velocities, or framelets which are located on the boundary but are moving 
tangentially to the boundary. 

To assess the time complexity of the TDPSF algorithm, we recall remark [^T^ 
which says that when Lint ^ w, the cost of computing Pout (by Algorithm 
I2.2p is only 0{L^^^k^^^log{Lintkma.x))- Thus, in this case, we find that the 
added complexity of TDPSF propagation over the Split Step algorithm is of 
lower order than the Split Step itself. However the constant is significantly 
larger in our experiments, so on a small grid the TDPSF propagator is slower 
than FFT/Split Step propagation. 

We discuss some possible improvements on this algorithm in section [5] 

3 Numerical Examples 

In this section we discuss the results of our numerical tests. 

The TDPSF algorithm is built in the program Kitty. Kitty is written in 
Python, with C extensions, calling the external libraries FFTW [12], Numarray 
and Matplotlib. The external programs gnuplot, ImageMagick and gifsicle were 
also used for making graphs/movies. 

Kitty is licensed under the GPL. It is very much a work in progress, and has 
little documentation and minimal user interface (an end-user version is currently 
in progress). Various test cases, spanning many types of parameters, are also 
available for download from the author's webpage, jhttp: / /math.rutgers.edu / "stucchio 

3.1 Simple Tests: Free Schrodinger Equation 

The standard method for testing absorbing boundaries is simply to throw coher- 
ent states (which are well localized in frequency) at the boundary and compute 
the difference between the approximate result and the true solution. This is a 
useful test, although it by no means completely characterizes the errors, as we 
discuss in Section [3T2l 

Our implementation is as follows. The free Schrodinger equation (i.e. 
g(t,x,^{x,t)) = 0) was solved on [-102.4,102.4] with TDPSF boundaries on 
the regions [-102.4,-88] and [88,102.4]. The lattice spacing was Sx = 0.1 in 
space, 5 x 10""* in time, and the TDPSF filtering interval was Tgtep = 2 x 10^''. 
The initial condition was taken to be ip{x, 0) = e~^^/^e™^ with v ranging from 
1 to 25. The simulation was run from t — out to Tmax — 3 • 51.2/'i;. 

After the solution was given sufhcient time to exit the computational domain, 
the error in the region [—88, 88] was measured (comparing with the exact result 
(|2.9p ). The error as a function of velocity is graphed in Figure [T] The TDPSF 
was used with cr = 1, 2, 4. The results are comparable to the complex absorbing 
potential V{x) = —25ie~'^^~^^^^'> which is also shown in Figured] The width 
of the complex potential was chosen so that it's spatial extent is comparable to 
the width of the TDPSF used. Altering the height of the potential changed only 
position of the dip in Figured] and achieved little other benefit. 
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5 10 15 20 

Pulse Velocity 

Figure 1: A graph of error vs the velocity of an outgoing pulse. The boxes, 
diamonds and "+" signs indicate TDPSF error reports. 

In step 2c of Algorithm 12.41 the TDPSF checks whether it is making an 
error. The above simulations were run with tolerance 10~^, and the simulations 
for which the TDPSF reported an error are indicated in Figure [T] Each time 
the TDPSF failed (and 4 times when it was successful), it notified the user of 
the error. The absorbing potential was unaware of it's own failures. 

This particular example demonstrates no major advantage of the TDPSF 
over the absorbing potential apart from the awareness of errors when the occur. 
The TDPSF works better for high frequencies, but not low. The advantage of 
the TDPSF is not that it succesfuUy dissipates outgoing waves, but that it does 
not dissipate incoming waves. An example of this will be demonstrated in the 
next section. 

Remark 3.1 Using the programs from [15J (available at'http:/ /www. math. unm.edu/~hagstrom/ 1, 
we ran a similar set of simulations to the ones above using the PML for bound- 
ary absorption. Since FDTD methods rather than spectral methods were used, 
we took 5x = 6.7 x 10"^, 5t = 2.4 x 10"^, and a PML with a thickness of 500 
lattice pointfl The error was of order 10"^ — 10"'^, probably due to the use of 
FDTD instead of spectral methods for solving the interior problem. 



^The computational region was [—20, 20] with 30, 000 lattice points. Solving on [—100, 100] 
would require require 5 X 30,000 = 150,000 lattice points (the TDPSF/spectral solver uses 
512 lattice points for the region [—20,20]). 
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3.2 Harder Tests: Long Range Potentials 

The main problem associated with an absorbing potential or PML is that not 
all waves located near the boundary are outgoing. The problem is that some 
waves are incoming, and should not be dissipated. 

Consider the following linear Schrodinger equation (with {x,t) G 



idti^{x,t) 



0.05|f|^ + 1 



V'(x,t) (3.1) 



^(X,0) = g^7.2g-|^|^/20^g.4xlg-|^|^/20 

The initial condition consists of two coherent states of equal mass, one with 
velocity 4 and one with velocity 7. The notable fact about this particular 
potential is that the fast Gaussian has enough kinetic energy to (mostly) escape 
from the binding potential, while the slow Gaussian does not. The slow Gaussian 
moves toward the boundary, turns around and returns. 

The problem with the absorbing potential approach is that it does not dis- 
tinguish between incoming and outgoing waves. It dissipates everything on the 
boundary including the waves that should have returned. We believe this should 
also occur for the PML (at least the PML found in [T5]). 

We ran three simulations of (|3.ip . The first was performed using the TDPSF 
with cr = 2.0, xo = 0.8, fco = 27r/12.8 = 0.491 and T^top = 0.1. The region of 
computation was [—25.6, 25.6]^, with lattice spacing 5x — 0.2 in space and 
5t = 0.025 in time. The region of interest was [—10, 10]^. Note that even if the 
region of interest were larger (as in Section [3?T|) . the width of the TDPSF region 
would not gro\\{^. An identical simulation was performed (on the same region) 
with an absorbing potential 

Fi(x) = -20ie-(^~i±25.6)V36 _ 20je-(--2±25.6)V36 

The third was solved with periodic boundary conditions on [—204.8, 204.8]^ 
(all other parameters were the same). This boundary is sufficiently distant so 
that the outgoing waves cannot return to the origin for a time 408.6/7.0 ~ 58. 
Letting ipt{x,t) be the solution with the TDPSF boundary, ■0a(a;,O be the 
solution with the absorbing boundary and ijjd{x,t) be the solution with the 
distant boundary, we measured the relative error: 

^ llV't,a(a;,t) -V'<i(a:,i)llL2([_ioio]2) 



\\M^)\ 



Figure [5] plots the error. In total, the error for the absorbing potential 
reaches 7.0% hy t = 58, while the TDPSF makes an error of 5.0%. The reason 
the error is so large for the TDPSF is that there is a concentration of mass near 



^We could have performed this test on the box [-204.8, 204. 8]^, with the TDPSF region 
still having width 15.6. But comparisons with an accurate simulation on a larger box would 
be impossible. 
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fc = near the boundary. This poses a problem for the absorbing potential, as 
Figure [2] shows. 

However, the problem with the absorbing potential is even more serious 
than it appears at first glance. The true solution can be written as 

ip{x,t) = J2j s~''^^*i^j(t>j{x) + fpd{x,t), with 4>j{x) the bound states and ipd{x^t) 
the dispersive part of 'ip{x,t). According to our simulations, by approximately 
t = 30, ipdix^t) appears to have dispersed. Therefore, after this point, all that 
should remain are bound states. 

This implies that M{t) = \\'^P{x,t)\\^2(^_iQ iq^) should remain approximately 
constant after t = 30. This behavior can be seen in the distant boundaries sim- 
ulation as well as the TDPSF simulation. This does not occur for the absorbing 
potential. Figure [3] graphs M{t) for t G [0,60] (in the distant boundaries case) 
and t £ [0, 120] otherwise. The absorbing potential appears to be dissipating, 
although we know this is impossible. For larger times, even the qualitative 
behavior of the absorbing potential simulation is wrong. 

The reason the TDPSF performs better than the complex potential is that it 
distinguishes outgoing waves from incoming waves. The TDPSF only removes 
waves which sit on the boundary and are also clearly outgoing. Low frequency 
waves on the boundary, which return to the computation region are incorrectly 
dissipated by the absorbing potential, but (correctly) ignored by the TDPSF. 
This is why the TDPSF gives the correct long time behavior. 

This problem is more dangerous than this example suggests. Although in 
this example we can determine that the dissipation is artificial, we cannot always 
do so. Similar problems occur often in atomic physics, but with exponentially 
decaying resonances as well as bound states (corresponding to complex Ej). 
In this case, one cannot distinguish the dissipation caused by the absorbing 
potential from that caused by the dynamics, and the measured decay rate will 
be incorrect. Since the object of such a simulation is often the measurement of 
the decay rate, this is a serious problem. 

3.3 Soliton Filtering 

Consider the nonlinearity, g{t, x,tp{x,t)) — — \4>{x,t)\^. It is desirable to con- 
struct a numerical algorithm which filters outgoing solitons as well as free waves. 
Although the TDPSF was not designed to filter solitons, it turns out to work 
well for solitons which are not moving slowly. This is however a coincidence, 
and cannot be expected to hold in general. A referee pointed out that it will 
fail for the KdV equation. 

The reason for this is that an outgoing soliton with sufficiently high veloc- 
ity is localized in phase space on outgoing waves. Consider a simple soliton, 
4>{x,v,t) = 2~-^/^e'("'^"'"(-^~'" cosh((a:: — vt)). A simple calculation shows that 
for |/c — w| 3> 1, (j){k,v,t) ^ e^l'^^"!. Thus, for k ^ fcmin, this shows that the 
framelet coefficients of (f){x,v,t) which are moving too slowly to resolve have 
exponentially small mass. This shows that under the free flow this soliton is 
strictly outgoing. 
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Figure 2: A graph of the error in 10, 10]^) for the absorbing potential and 

TDPSF up to i = 58. 



Local Norm Measurement 



18 
17 



M{t) 



' TDPSF 
Distant Boundaries 
Absorbing Potential 




Figure 3: A graph of M{t) = \\^j{x,t)\\ j^2([_i(f xo]2) ^- The distant boundary 
simulation is invalid at time t = 58, due to the return of outgoing pulse. 
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Figure 4: A plot of the error (defined as in ()3.2p ) as a function of velocity. Note 
the exponential improvement in accuracy with velocity. 

The soliton is also leaving the box under the full &owU{t). Although e*'^^/^)'^* 
and U{t) move the soliton very differently (one dispersively, one coherently), 
they both move it out of the box and in nearly the same direction. For this 
reason, we expect the TDPSF is to filter soliton solutions correctly. 

We ran numerical tests to demonstrate this as follows. We solved (|l.ip 
with g{t,x,tp{x,t)) = — \ip{x,t)\'^ on the region [—25.6,25.6] with Sx — 0.05, 
6t = 0.002/z; and Tstop = 0.008/u (the timestep's are scaled with the velocity to 
speed up the simulations). In this simulation, Lint = 12.0 and w = 13.6. The 
initial condition was taken to be iIj{x,0) = 2^^/^e™^/ cosh(x) for v — 1..15. 

The TDPSF was used with Tstcp = 0.08/u, xq = 0.20, fco = 27r/3.2, and 
(7 = 1,2, 3. We measured the following quantity: 

, llV'(a;,t) - V'ex(a;,t)||i2([_i2,i2]) ,^ ^, 

E{v) = sup II , f (3.2) 

t<200/t) \\y^ex[x,V}\\j^2^g^^ 

The function ipex{x,t) is the exact solution. The result of this experiment is 
plotted in Figure [3T3l The time 200/ti was chosen since it is more than enough 
time for errors to return to the region [—10, 10]. We believe the error floor near 
10"^*^ visible in Figure [3731 is due to truncation errors in the calculation of the 
WFT, and errors due to time-stepping in solving the interior problem. 

Remark 3.2 The paper [26] proposes an alternative method of absorbing bomid- 
aries (namely the paradifferential strategy), based on a novel method of ap- 
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proximating the Dirichlet-to-Neumann operator. A similar numerical test was 
performed for those boundary conditions. For a soliton at velocity 15, Szeftel 
obtained £'(15) = 0.08 at best. For comparison, we obtain E{15) = 1.69 x 10~^° 
for cr = 1 and E{v) = 1.40 x 10"!° for cr = 3. 

Our tests are not directly comparable. Our region of interest was [—12, 12] 
with TDPSF region on [-25.6,-12] and [12,25.6], as opposed to Szeftel who 
used [—5, 5] (although we used 256 spatial lattice points rather than 200, as 
used in [55] )■ We used FFT/Split Step propagation for the interior problem, as 
opposed to the finite differences of [55]. Additionally, storage of the history on 
the boundary was required in [26j . unlike the TDPSF. 

It is quite surprising that the TDPSF beats the Dirichlet-to-Neumann bound- 
ary by such a large margin because |26| takes the nonlinearity into account while 
the TDPSF assumes the nonlinearity is zero on the boundary. 

4 Comparison to Other Methods 
4.1 Dirichlet-to-Neumann Boundaries 

The Dirichlet-to-Neumann map was originally constructed by Engquist and Ma- 
jda in TT,]!] (see also [l][2])- Their guiding principle was that near the boundary, 
the physical optics approximation to wave flow is sufficiently accurate to filter 
off the outgoing waves. The TDPSF is a direct analogue of this - the Gaussian 
framelet elements behave (under the free flow) like classically free particles. We 
use a different method to filter, but the guiding principle is the same. 

Modern approaches attempt to construct the exact solution on the boundary 
and then impose it as a boundary condition. In principle, this is the best 
possible approach, although in practice it is difficult to implement for dispersive 
equations. We briefly mention two major approaches that we are aware of, and 
remark that only one [25] even attempts to deal with nonlinear equations. 

An additional problem is that this approach precludes the use of spectral 
methods (e.g. Algorithm [2T3|) to solve the interior problem, reducing accuracy of 
the solution on the interior for dispersive equations. The FFT naturally imposes 
periodic boundaries rather than Dirichlet-to-Neumann, thus requiring the use 
of FDTD or other local methods. 

4.1.1 Current Constructions 

To deal with the free Schrodinger equation (no nonlinearity or potential) , Lubich 
and Schadle [151 US] constructed an approximation to the exact integral kernel 
of the Dirichlet-Neumann operator based on a piecewise exponential (in time) 
approximation. This approach appears to work nicely for the free Schrodinger 
equation, although it is uncertain that it could be applied to the full Dirichlet- 
to-Neumann operator of a nonlinear equation or long range problem. 

We are aware of only one fully nonlinear Dirichlet-to-Neumann operator, 
constructed by J. Szeftel in [25]. Szeftel uses a modified version of the parad- 
ifferential calculus (see references in [3S]) in 1 space dimension to deal with 
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smooth nonlinearities and potentials. He proves local well posedness of of 
the boundary conditions assuming regularity of the initial data. 

However, extensions to appear highly nontrivial. The assumptions are 
significantly stronger than ours, and there are no error bounds. The numer- 
ical experiments look promising and the results appear accurate for radiative 
problems (see also Remark l3.2p . 

4.2 Absorbing Potentials/ PML 

4.2.1 Absorbing Potentials 

Absorbing (complex) potentials, described in [T7], are the current "industry 
standard". One can add a dissipative term —ia{x)'ip{x,t) to the right side of 
(jl.ip and solve it on the region [— (ijnt + w), (iint + w)]^ . The function a{x) 
is a positive function supported in [—(Lint + w), (^int + w)]^ \ [— Lint, -^int]^- 
This has the effect of (partially) absorbing waves which have left the domain 
of interest, although it might create spurious reflections and dissipation. This 
approach is the mainstay of absorbing boundaries, being simple to program and 
compatible with spectral methods. 

The potential a{x) must be tuned to the given problem. Given fcmin, fcmax, 
one must select the height and width of the absorber so that it kills most of the 
wave between fcmin and fcmax- Waves with momentum lower than fcmin are mostly 
reflected, and waves with momentum higher than fcmax are mostly transmitted 
and wrap around the computational domain. 

Heuristic calculations and numerical experiments suggest that the absorber 
must have width proportional to Cfcmax ln(e) / fcmin, with C depending on the spe- 
cific shape of the potential. The TDPSF works on a layer of width C ln(e) /fcmin, 
which is smaller by a factor of fcmax- Spurious dissipation is an additional prob- 
lem with absorbing potentials, as illustrated in Section [5?^ 

4.2.2 Perfectly Matched Layers 

Perfectly Matched Layers (PML) are a variation on this approach, first proposed 
in [3] for Maxwell's equations. In |T5j they are designed and tested for the 1 
dimensional free Schrodinger equation, with reasonable results. 

In the special case of the Schrodinger equatior0, the PML consists of adding a 
term —ia{x) (1/2) Atp{x,t) to the right side of (jl.ip instead of merely —ia(x)?/'(a?, t) 
If a{x) is chosen carefully, one can completely eliminate reflections at the in- 
terface (the boundary of suppa(x)) for certain frequencies. It also has the 
property of dissipating high frequency waves more strongly than low frequency 
ones, thereby requiring a boundary layer of size only Cln(e)/fcmin. 

As discussed in remark [5TT1 numerical tests using the PML of [15] achieved 
errors of only 10^'^. The PML is likely to have the same problem as complex 

•^In general, one constructs the PML by exterior complex scaling, but for the Schrodinger 
equation it takes this simple form. 
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absorbing potentials with spurious dissipation. Some PML methods are un- 
stable [18] . The PML method for the Schrodinger equation is still very much 
undeveloped, so a more detailed comparison is difficult to make. 

5 Conclusion 

We have described in this work a new method, the Time Dependent Phase 
Space Filter, of filtering outgoing waves for the nonlinear Schrodinger equation. 
Unlike absorbing potentials and (most likely) PML methods, the TDPSF filters 
only those regions of phase space containing outgoing waves. Dissipative terms 
on the boundary filter all waves, including waves which should have returned to 
the computational region from the boundary (see Section 13. 2[) . 

Our method is easier to construct than the Dirichlet-Neumann map, be- 
ing based on standard techniques of signal processing rather than complicated 
pseudo/para-differential calculus. It is local in time, whereas the Dirichlet- 
Neumann map is history dependent. The TDPSF is also compatible with spec- 
tral methods, unlike the Dirichlet-to-Neumann map. 

Unlike all other methods we are aware of, the TDPSF fails gracefully, in the 
sense that when it cannot filter the outgoing waves it notifies the user rather 
than reporting an incorrect calculation as correct. 

The main problem with the method is the problem of filtering low frequency 
outgoing waves, a problem shared with all other methods. We are currently 
working on an extension of the TDPSF [23] which will reduce the cost of filtering 
from Cln(e)/fc,jii,i to Cln(e) In(fciiiin). 
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NSF grant DMSOl-00490. 
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